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ABSTRACT 

Mechanically induced viscoelastic dissipation is difficult to compute. When the constitutive 
model is defined by history integrals, the formula for dissipation is a double convolution integral. 
Since double convolution integrals are difficult to approximate, coupled thermo-mechanical 
analyses of highly viscous rubber-like materials cannot be made with most commercial finite 
element software. In this study, we present a method to approximate the dissipation for history 
integral constitutive models that represent Maxwell-like materials without approximating the 
double convolution integral. The method requires that the total stress can be separated into 
elastic and viscous components, and that the relaxation form of the constitutive law is defined 
with a Prony series. 

Many commercial codes employ history integral formulations to compute viscous stresses. 
The numerical approach taken to approximate the history integral often involves an interpolation 
of the history integral's kernel across a time step. Finite difference equations are then derived for 
the evolution of the viscous stresses in time. When the material is modeled with a Prony series, 
the form of these finite difference equations is similar to the form of the finite difference 
equations for a Maxwell solid. This analogy implies that the value of a Maxwell solid's 
dissipation function may be close to the value of the dissipation function of a history integral 
constitutive model. Since the dissipation rate in a Maxwell solid can be easily computed from 
knowledge of its viscous stress and Prony constants (spring-dashpot constants), we investigated 
employing the Maxwell solid's dissipation function to couple thermal and large strain mechanical 
finite element models of rubber. Numerical data is provided to support this analogy and to help 
understand its limitations. Rubber cylinders with imbedded steel disks and with an imbedded 
steel ball are dynamically loaded, and the nonuniform heating within the cylinders is computed. 


INTRODUCTION 

Rubber is employed to carry large loads in tires, gaskets, and tank track pads. It is also used 
to provide damping and system stability in complex mechanical systems such as helicopter 
rotors. In these applications, the rubber is typically stiffened by the addition of carbon black. 
The filled rubber tends to be a poor conductor of heat, yet it also exhibits very large hysteretic 
energy loss during cyclic loading. Also, the mechanical properties of rubber are strongly 
dependent on temperature. Faced with the above issues, designers interested in modeling the 
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detailed response of complex- shaped viscoelastic rubber components need to be able to compute 
the coupled thermo-mechanical behavior of rubber. 

It is computationally difficult to compute dissipation in viscoelastic materials modeled using 
history integral methods. The dissipation function is expressed as a double convolution integral. 
Commercial finite element codes that do not compute dissipation for materials with memory lack 
the ability to perform coupled thermo-mechanical analyses of rubber components. 

In this paper, an approximate dissipation function is introduced for viscoelastic materials 
whose relaxation moduli are modeled with a Prony series. The dissipation function is 
numerically evaluated for a range of strain rates considered as fast to slow with respect to the 
Prony series decay rates. The dissipation function is employed to couple the large deformation 
rubber viscoelasticity and heat transfer finite element algorithms. The resulting thermo- 
mechanical algorithm is used to simulate the hysteretic heating of dynamically loaded rubber 
cylinders. Again, the purpose of the work is to establish and test a computational tool for 
analyzing hysteretic heating in rubber components. 

Overviews of large strain rubber viscoelasticity and heat transfer are included to facilitate the 
description of the thermo-mechanical coupling performed in this study. Only moderate 
temperature changes were simulated in this study, so time-temperature superposition is not 
discussed in this effort. 


DISSIPATION IN MATERIALS WITH MEMORY 


The thermodynamics of materials with memory was formulated and investigated in detail by 
the early 1970s. We reference a few papers here to provide an overview on the computation of 
dissipation in materials with memory. Stress wave propagation in viscoelastic solids was 
investigated in 1961 by Hunter [1]. In his work, Hunter included a discussion of the 
formulations for thermo-mechanical coupling. A technical presentation of the viscoelastic 
dissipation in its double convolution integral form was also included. He analytically 
investigated the coupling effects by studying the case of steady state harmonic strains. Complex 
modulus methods were then employed to produce a dissipation function. 

Several papers by Oden and Aguirre Ramirez [2], Oden and Armstrong [3], and a monograph 
by Oden [4] present finite element algorithms for thermo-mechanical analyses of materials with 
memory. In these references, Oden, et al. present the finite element formulation for coupling the 
equation of motion and the heat conduction equation for materials with memory. To numerically 
accomplish the coupling, Oden, et al. approximated the viscoelastic dissipation (the double 
convolution integral) for the case when the material relaxation is represented with a Prony series. 
They employed finite difference methods for the strain rates and numerical integration for the 
double convolution integrals. 

Using the notation of Hunter [1], the stress in a linear one-dimensional solid can be described 
as follows 


o 


(<)=£ e(t)~ 



( 1 ) 


where o is the stress, e is the strain, t is the current time, and E is the material’s elastic 
modulus. The function 0 is of the form 
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in which g(z)> 0. The viscoelastic rate of energy dissipation, r(/), for this linear solid (see 
Hunter [1], Oden [4], and Christensen [5]) is expressed with the following double convolution 
integral 


r(t)=E J J ) jL z)(2t - si - s 2 )ds x ds 2 > 0 

J J dsi 0 S 7 at 

s 2 =—°°s l =—°° 

Generally, this dissipation function is difficult to evaluate numerically [2, 3, 4]. 


(3) 


APPROXIMATION OF DISSIPATION FOR MAXWELL-LIKE MATERIALS 


To avoid the computation described in Equation (3) we exploit the fact that when the 
relaxation modulus in the history integral constitutive model is in a decaying Prony series form 
(that is, the function Q in Equation (2) is described with decaying exponentials) then the history 
integral can be approximated with a finite difference equation that is similar in form to the finite 
difference equation for a Maxwell solid. We explore the use of the Maxwell solid's dissipation 
function (which is easily computed from stresses) for approximating mechanically induced 
heating in viscoelastic solids defined by history integrals. 


Dissipation Function for a Maxwell Solid. In a Maxwell solid (see Figure 1) with a spring- 
dashpot leg having a spring modulus of E v , a spring strain of £ 1 , a dashpot with viscosity of rj , 
and a total dashpot leg strain of £ , the rate of energy dissipation r(t ) is given as follows. 


dissipated work 
unit time 



(4) 
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de de 


or 


r(t) = 


dt dt 


> 0 


n 




with E v £ v = 77 




d\s - 
dt 


V 


in the spring-dashpot leg. Equation (5) becomes 


(5) 
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where U v (t) is the energy in the spring on the spring-dashpot leg (the viscous leg) and T = 

is the relaxation time. When the rate of dissipation is measured per unit volume, the term U v (t) 

is the corresponding energy density. Then, values of U v (t) can be approximated by the current 

values of the viscous stress, <7* (/). If the spring-dashpot leg has a time constant given by l , the 
resulting rate of energy dissipation per unit volume is given by 




If we assume the Maxwell model has a rubber spring-dashpot leg, and if we compute stretch 
ratios from knowledge of the viscous stresses, then we can approximate the linear solid’s energy 
density U 1 (/) with a rubber energy density function and use it in Equation (6). 


Finite Difference Equations for Viscous Stress Components. A simplified one-dimensional 
version of the finite difference equations (across a time interval At = t n+ \ - t n ) for the evolution 

of the viscous stress, <7 v (t), in a history integral model employing Prony series material 

constants is given by [6] 



where (jjj is a stress increment computed at time l n by using the relative deformation gradient 
between configurations at times t n _\ and t n , and using the instantaneous modulus, Eq , of the 

material. The analogous finite difference equation for a Maxwell link (obtained by employing 
the trapezoidal method) is given by 


4 



-At 

>g + e T o v n 


(9) 


®n + 1 ^ 


^0 0 
^n+1 


f At ^ 

1 + — 
2t 


where ((7° +1 -(7°) is a stress increment computed using the total strain (i.e., current 
configuration's strain) and using the instantaneous modulus of the material. 


For a time interval in which the stress increment terms, {• • , in Equations (8) and (9) are 

approximately equal, or in which the stress increment terms have small values relative to the 
-At 

stress decay term, e 1 <j]\ , the Maxwell model and the history integral model will have stress 

histories that are approximately equal. In this case, the dissipation functions for the two models 
will yield approximately equal values for that time interval. 


Procedure for a Large Strain History Integral Viscoelasticity Code. In this section we 
provide an outline of a variational formulation for rubber elasticity, a history integral version of 
rubber viscoelasticity [6], and a procedure for approximating the viscoelastic dissipation. 
Notational definitions that are not included here are given in the Appendix. The coordinates of a 
material point in the reference configuration are indicated by X = {A, , A 2 , A 3 } and in the 

deformed configuration by x = {xj , x 2 , x, }. The vector function x = x(X, /) defines the mapping 

of points from the reference configuration into the deformed configuration. The Lagrangian 
method is employed to describe the deformed solid's energy. The virtual work statement, 
excluding inertial effects, for a deformed solid of volume V and surface area S is 

SW I = J<r : <SD V dV = J Sv T t dS + J Sv T f dV (10) 

v s v 


where 6\ is a virtual displacement vector, S\V I is the internal energy due to the virtual 
displacement, t is the traction stress prescribed over S , f is the body force vector acting within 
the volume V , and V 0 is the volume in the undeformed state. When the solid is rubber, the 
internal energy in Equation (10) is expressed as follows. 
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(11) 


where U is the hyperelastic strain energy density function defined in this effort as 


U = ^-(h-3) + U(J- 1 f 


( 12 ) 
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and where G is the shear modulus, and K is the bulk modulus. 

Quasi-static deformations of a viscoelastic rubber solid are often modeled by approximating 
the material constants with a Prony series. That is, the shear and bulk moduli are expressed in 
the form 



N —■ 



N — 

G(f)=G 0 

^0 

[XI 

+ 

8 

&0 

and 

K(l)=K 0 

XI 

+ 

8 


i — 1 



. ' =1 J 


(13) 


where Gq = G(t = 0) and K 0 = K(t = 0) are the instantaneous shear and bulk moduli, and the 

D hf 

Xj ’s denote relaxation time constants. The deviatoric, a , and hydrostatic, a , stresses are 
determined from history integrals as follows. 



and 

N , 

dt (15) 

i = 1 1 

where F r (t - E , ) is the deformation gradient of the deformed shape at time t-E, relative to the 

deformed shape at time t . As mentioned above, approximations are made when Equations (14) 
and (15) are integrated so that the entire history does not have to be convoluted. The 
approximations depend on the material constants being of Prony series form. Throughout the 
remainder of this effort only one term in the Prony series is employed, and the Prony series 
subscripts are dropped. 

The following procedure was employed (at each integration point in each finite element) in 
this study to compute the approximate viscoelastic rate of dissipation (Equation (6)): 

Procedure. 

5x 

A. ) Kinematics . Obtain the deformation gradient, F = — , from the finite element (FE) code. 

OX 

Compute the volume change measure, J = det(F) , compute the deformation gradient with the 
volume change eliminated, F = 7 1/3 F, compute the deviatoric stretch matrix (the left Cauchy- 
Green strain tensor), B = F ■ F T , and find the first strain invariant, /, = trace B . 

B. ) Stresses . Obtain the Cauchy stresses, < 7 , from the FE code. Compute the hydrostatic stress, 
p = —— I: <r, compute the deviatoric stress S = <7 + /; I , compute the elastic deviatoric stress, 
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. Note, in this effort, we employed an energy density 


S E = -DEV 
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dU y dU 
d/ i dl 2 


|b-^b 2 

3/-, 


function U that was not a function of I 2 . Then, since /, = trace B , our calculation becomes 

f T \ 

for the normal stresses and 


S E =jDEv[g e b]=j 


G e B - — trace(G E B ) 


= -g e 

J E 
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3 


V 
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S' = — G /. B for the shear stresses where G E = Gq( 1 - g). The deviatoric viscoelastic stress is 


/ 


computed next using S v = S - S E . Compute the principal viscoelastic stresses, , S v p2 , S p 3 J, 

based on S' . In the case when a rubber energy density function will be used, compute the 
principal stretches for the viscoelastic part using the following equations [7]: 


r* V o V 

b p\ ~Jp2 

V V 

Jp2-Jp3 

nV ri V 

*p 3 ~ ^ pi 



and the constraint equation X- ■ X v • X,, = 1 • To solve for the stretches, we converted these 
equations into a cubic equation of the form 


x 3 + r ■ x° + s ■ x + 1 = 0 

where x-X v , r = (xj v ~ X 2v )= - 5^ 3 )/ G v , 5 = -(/L 2 V - X 2v )(X> V - A 2 V )= 

- (s' ^ | - Sp 2 fs p2 - S p2l )/ G 2 , and t = -1.0 . Following procedures for solving cubic 

equations, the equation can be reduced to 


y 3 + p-y + q = 0 


2 2 o r ■ s 

where y - x + r/3, p = s-r / 3 and q = --r — — + 1 . By solving the above cubic equation 

and substituting y into the previous equation, the principal stretches (A] v , X 2v , Xv ) can be 
obtained. 


C.) Dissipation . In the case when the Maxwell solid is assumed to have a rubber spring-dashpot 

2 2 2 

leg, compute I\ v = \ v + Xv + Xv • Then compute the viscoelastic energy density with 


,, G,, ,, 2 U v (f) 

U = — - (/ 1 v , - 3). The rate of energy dissipation is then approximated with r(t)= . In 

2 T 
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the case when the Maxwell solid is assumed to have a linear spring-dashpot leg, use the viscous 


deviatoric stresses to compute the rate of dissipation as 



S V :S V 
4 Gy T ' 


Accuracy of the Dissipation Function for Rubber. The approximations to r(t ) described 
above are evaluated by solving benchmark problems. We select several uniform axisymmetric, 
plane stress, and plane strain problems (see Figure 2) and determine the dissipation by 
calculating the work done during closed load-displacement cycles. Three measures of the 
dissipation are calculated. The cyclic work done by the finite element code’s external forces, 
and the total dissipated work found by integrating Equations (6) & (7) over the solid’s volume 
and over time. The time frames and strain rates in these benchmark examples are selected 
relative to the Prony series relaxation time 1 as follows: 

Case 1 . A generic description of the first cyclic loop employed is given in Figure 3. The loop 
represents a rapid step-strain loading (0,q), followed by a relaxation (q ,t 2 ), and then followed 
by a reduction in strain and load until the loop was effectively closed (/ 2 ^ ^ 3 ) • 


Case 2 . A generic description of the second cyclic loop employed is given in Figure 4. The loop 
represents a rapid step-strain loading (0, q ), followed by a relaxation (q ,/ 2 ), followed by a rapid 
step- strain unloading to zero load (/ 2 ’ ^ 3 ) ’ an d ending by relaxing until the loop was effectively 
closed ((3, 4). 

Case 3 . Dynamic tension-tension and compression-compression enforced sinusoidal 

displacements were applied at three frequencies. The frequencies were 

co =0.1Hz, 1.0Hz, and 10.0Hz (slow to fast relative to the material’s viscous relaxation time 
t - 1.0 sec ), see Figures 5 to 7. 

The ABAQUS [6] finite element code was employed to determine the work done in closed 
cycles for each Case listed above. The four node rectangular CAX4HT axisymmetric, CPE4HT 
plane strain, and CPS4T plane stress elements were used. The elastic and viscous material 
constants used in this accuracy study are Gq=1.0 MPa, g = 0.5, t= 1.0, and 

Kq = K rX} =3000. MPa. The thermal material constants are listed in Table 1. The finite 

element code’s reaction forces were used to compute the total internal dissipation for each of the 
cycles. The results are referred to here as J F T du . The Procedure listed above for 

approximating the dissipation was employed to compute the total dissipation in a cycle. The 
approximate dissipation results are referred to here as f r(t)dt . The Error is defined as 


Error = 


J F T du - jr(t)dt 



(16) 
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The error calculations for axisymmetric loading produce data nearly identical with the error 
data obtained for plane stress. The errors obtained in plane strain are larger than those obtained 
in plane stress when they are compared relative to the maximum principal strain. 

In the case of axisymmetric loading, the accuracy results are listed in Tables 2, 3, and 4. The 
approximation to the dissipation has an error less than 2.0% for small strain deformations and the 
error grows large as the strain amplitudes are increased. At maximum tensile strain of 0.10 in 
Cases 1 and 2, the errors in the approximation to the dissipation are approaching 10%. 

At large strains, the errors are large, and corrections are needed if this approach is used. It 
appears possible to calibrate the error at large strain. For example, Figures 8 and 9 show the 
error as a function of the maximum principle strain and as a function of the first strain invariant 
for both plane stress and plane strain. A look-up-table could be constructed based on strain 
invariants or stretch ratios to provide correction factors at large strains. 

In these examples, approximating the dissipation with the formulas for energy in a rubber 
spring-dashpot Maxwell leg yielded smaller errors than when the energy formula for a linear 
solid-dashpot Maxwell leg was used. In the nonuniform thermo-mechanical heating examples 
presented below we did not adjust the approximation to the dissipation with a look- up-table. 


THERMO-MECHANICAL HEATING OF RUBBER SOLIDS 

Heat Transfer. The variational statement of the energy balance equation for heat transfer, 
together with Fourier's law, for a deformed solid of volume V , and surface area S is given by 

f p^-50dV + k ^ dV = f S0rdV+ \dOqdS (17) 

J dt J dx dx J J 

V V VS 

where 6 is the temperature, U e is the internal thermal energy, p is the density, q is the heat 
flux per unit area, r is the heat supplied per unit volume (the dissipation described above), k is 
a conductivity matrix, and 56 is a virtual temperature field satisfying the essential boundary 
conditions. Equation (17) is used to build the transient heat transfer finite element equations for 

the rubber solid. The term j 50 r dV couples Equation (17) to Equation (10). 

v 

Finite element discretization of the mechanical and thermal variational statements, Equations 
(10) and (17), results in systems of time dependent differential equations which approximate the 
mechanical and thermal response of the solid. Below, we computationally investigate heat 
buildup in rubber solids by employing the approximate calculation of viscoelastic energy 
dissipation described in the Procedure above. 

Finite Element Employed. We present the results for the heating of several 2D uniform and 
non-uniform axisymmetric models. The eight-node hybrid axisymmetric CAX8RHT element 
was used, see Figure 10. The interpolation in the element is biquadratic for the displacement, 
and bilinear for the temperature and pressure. Reduced integration is also employed. The nodal 
variables employed in the thermal element are the temperatures, 6 , at the four corner nodes. 
The stress element has radial, u r , and axial, u, , displacements at the eight nodes and hydrostatic 
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pressure, p t , variables at each of the four integration points. The computed heating rates, r(/), 
are applied in the heat transfer algorithm at the four integration points. 

Five cylinders were analyzed for hysteretic heating. Only moderate temperature changes 
were being simulated so the elastic material constants were not treated as temperature dependent 
and time-temperature superposition was not employed. 

Cylinder Dimensions. There were five different cylinder models. The first four models 
consisted of two groups of two cylinders each. One group consisted of uniform cylinders and the 
other group had steel disks at their centers, see Figure 11. The first four cylinders each had a 
radius of R - 0.0282 m . There were two cylinder heights in each group. The heights were 
H = 0.05 m, and 0.0125 m respectively. The cylinders were compressed between steel fixtures. 
The model simulates the case when a lubricant maintains the fixture-rubber interface as 
frictionless. The internal steel disks were completely attached (bonded) to the rubber. The 
height and radius of the disks were 0.0025 m , and 0.0141 m respectively. The fifth model was a 
tall cylinder of height 0.10m. This tall cylinder had a steel ball imbedded at its center. The 
radius of the rubber cylinder was 0.05m and the radius of the steel ball was R=0.025m. The 
boundary conditions were the same as those used for the imbedded disk models. 

Rubber and Steel Properties. The elastic constants used for the steel are; Young’s modulus 
E- 206.8 G Pa, and Poisson’s Ratio v =0.3. The rubber energy density was modeled as a 
Neo-Hookean [7] solid, Equation (12). The viscous behavior was described with one Prony 
series term, Equation (13). The material constants employed are representative of filled rubber 

/n 7 7 

[8, 9], They are; — — = 1.155 MPa , = — = 1.000 * IQ 3 MPa , and 

2 2 2 

g = 0.3, k = 0.0 , z = 0.1 s . The thermal material constants for both the steel and the rubber are 
listed in Table 1. The film heat transfer coefficients for the rubber- air and rubber-steel interfaces 
are h A = 5.44284/ /(° C m 2 s) and h s = 20934 / /(° C m 2 s), respectively. 

Enforced Displacement. The deformation of the tall (H - 0.05 m ) uniform cylinder subjected 
to an axial displacement is shown in Figure 12. Since the boundary conditions are symmetric, 
only the half height (r, z)-quadrant is shown. The mesh refinement near the top and outer edges 
is to accommodate the heat transfer gradients at those locations. Each cylinder was subjected to 
an axial enforced displacement that produced moderately large strain hysteresis. The half height 
enforced axial displacement for the tall cylinder is described as follows. The top end of the 
cylinder was ramped to a displacement of - 0.0045 m in 0.05 s . The top end of the cylinder 
was then forced to follow the prescribed axial displacement given by 

u z (t) = -0.0045 - 0.0030 Sin(ln (6.5 )t) m (15) 

The half height enforced axial displacement for the short cylinder is described as follows. The 
top end of the cylinder was ramped to a displacement of - 0.001 125 m in 0.05 s . The top end of 
the cylinder was then forced to follow the prescribed axial displacement given by 

m,(z)= -0.001 125 -0.00075 Sin(2n(6.5)t) m (16) 
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where t is the time in seconds. The loading and the displacement of the top of the cylinder are 
shown in Figure 13 for a time interval of 1.0 s. As expected, the displacement curve 
demonstrates viscoelastic softening. A number of analyses were performed to test the nonlinear 
elastic, the viscoelastic, and the transient heat transfer finite element algorithms separately. 
Hand calculations and finite difference calculations verified that the algorithms functioned 
correctly. The thermal boundary conditions did not significantly affect the temperature fields 
computed. The following analyses were performed, employing the dissipation function 
described previously, to investigate the uniform and non-uniform heating of rubber cylinders 
undergoing large strain dynamic deformations. 

Uniform Cylinders. The cyclic deformations given by Equations (15) and (16) were applied to 
the tall and short uniform cylinders indicated by the dimensions given previously. The 
temperatures, at the center of the cylinders (r = 0, z = 0) as a function of time for the first 20 s of 
loading are shown in Figure 14. The frictionless rubber- fixture boundary condition allows the 
strains to be uniform in the cylinder. Rubber is a poor conductor of heat, and the fact that the 
cylinders heated nearly uniformly regardless of height was expected. 

Cylinders with Imbedded Disks. The models for the tall and short uniform cylinders described 
above were modified by imbedding steel disks at their centers. The models were cyclically 
loaded to simulate nonuniform viscoelastic heating in rubber solids when the deformations 
produce moderately large strains. The results obtained appear reasonable. Figure 15 shows the 
meshes on the reference and deformed shapes for the tall cylinder. The temperature distribution 
in the tall block after 20 s of dynamic loading is shown in Figure 16. Temperatures at points A, 
B, C, and D in Figure 16 are plotted as a function of time in Figure 17. The outer radial end of 
the internal disk (point C) is predicted to heat much faster than the rest of the cylinder. 

Similar results were obtained for the short cylinders. However, there were some convergence 
problems that were overcome by re-meshing near the top outer corner of the steel insert and 
employing lower order displacement interpolations (a four node displacement element was 
employed.) The temperature distribution in the short cylinder is shown in Figure 18. Again, 
Figure 19 shows that the region of moderately large strains, located near the outer radial end of 
the internal disk (point C), has the most rapid rise in temperature. 

Cylinder with Imbedded Ball. A tall rubber cylinder (H=0.10m) was modeled with an internal 
steel ball, see Figure 20. The properties and boundary conditions are the same as those 
employed for the models with the imbedded steel disks described previously. Note, this cylinder 
is much longer than the cylinders analyzed above and the overall deformation of this cylinder is 
less severe than the cases analyzed above. The top end of the cylinder was forced to follow a 
prescribed axial displacement given by Equation (15). The reference and deformed meshes are 
shown in Figure 21. The strains are moderately large. A temperature distribution profile after 20 
s is shown in Figure 22. The region above the ball near the surface of the cylinder is predicted to 
heat the fastest. The temperatures at these locations as a function of time are shown in Figure 23. 
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CONCLUDING REMARKS 


In dynamically loaded rubber structures predictions of both the strain and the temperature 
distributions are required to perform degradation studies. A procedure that couples a viscoelastic 
large deformation stress analysis with a heat transfer analysis to model nonuniform hysteretic 
heating in rubber was described. The ABAQUS finite element code was used to implement the 
procedure. 

Computations were made to demonstrate the thermo-mechanical heating of tall and short 
uniform rubber cylinders. The viscoelastic material properties employed are valid for 
moderately large strain deformations of rubber. When the time dependent strains in the cylinders 
were uniform, the heating computed was also uniform. The thermo-mechanical heating of tall 
and short rubber cylinders, containing internal steel disks and a steel ball were also computed. 
The internal steel components provided high strain gradients within the rubber cylinders and 
nonuniform hysteretic heating was observed. The analyses performed suggest that the coupling 
procedure should be considered for further development as a design tool for rubber degradation 
studies. 


REFERENCES 

1. Hunter, S. C., “Tentative equations for the propagation of stress, strain, and temperature 
fields in viscoelastic solids,” J. Mech. Phys. Solids, Vol. 9, pp. 39-51, 1961. 

2. Oden, J. T., and Aguirre Ramirez, G., “Formulation of general discrete models of 
thermomechanical behavior of materials with memory,” Int. J. Solids Structures, Vol. 5, pp. 
1077-1093, 1969. 

3. Oden, J. T., and Armstrong, W. H., “Analysis of nonlinear, dynamic coupled 
thermoviscoelasticity problems by the finite element method,” Computers and Structures, 
Vol. l,pp. 603-621, 1971. 

4. Oden, J. T., Finite Elements of Nonlinear Continua, McGraw-Hill Book Company, Library 
of Congress number 70-154237, 1972. 

5. Christensen, R. M., Theory of Viscoelasticity, An Introduction , Academic Press, Library 
of Congress Number 74-127681, 1971. 

6. Hibbitt, Karlsson & Sorensen Inc., ABAQUS Theory Manual, Version 5.8, 1080 Main St., 
Pawtucket, RI., 2000. 

7. Treloar, L. R. G., The Physics of Rubber Elasticity, Clarendon Press, 

ISBN 0 19 851355 0, 1975. 

8. Clark, S. K. and Dodge, R. N., “Heat Generation in Aircraft Tires Under Braked Rolling 
Conditions,” NASA Contractor Report 3629, 1982. 

9. Pitts, D. R. and L. E. Sissom, L. E., Heat Transfer, McGraw-Hill Book Co., Library of 
Congress number 77-20255, 1977. 


12 



APPENDIX 


The following notation is included to assist the reader with the description of the finite element 
algorithm for rubber viscoelasticity. 


Xf (i = 1,2,3) 
x { (i = 1,2,3) 
x = x(X,t) 

F = i* 
ax 


coordinates of a material point in the reference configuration, 
coordinates of a material point in the deformed configuration, 
vector mapping between the reference and deformed configurations. 

deformation gradient. 


7=det(F) 

_ 1 

F = /~3 F 

B = F F 7 


/>i((/7)-, r (BB)) 


determinate of F which measures volume change. 

deformation gradient scaled for volume change. 

left Cauchy Green strain tensor. 

first strain invariant (adjusted for volume). 

second strain invariant (adjusted for volume). 


<5u 

SL = 


dSu 

dx 


displacement, 
gradient of displacement. 


sd=^(sl+sl t 

SD y 

A : B 

Se vo1 =1 SD 


rate of deformation. 

rate of deformation computed using virtual displacement dv . 
scalar product of matrices A and B . 
volumetric strain rate. 


& = SD--Se vol I 
3 


P = 



: a 


deviatoric strain rate, 
pressure stress (hydrostatic). 
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TABLES 


Table 1. Thermal properties. 


Property 

Rubber 

Steel 

Conductivity, J / (°C m s) 

1 K = 

0.20934 

45.83379 

Density, kg / ni 

P = 

1000. 

7849. 

Specific heat, J /{kg °c) 

C = 

2093.4 

460. 

Expansion, ( 0 Cy 

a = 

80.*10“ 6 

12.* 10 6 


Table 2. Errors obtained predicting dissipation, Case 1, see Figure 3. 




2 U v (t) 

r(t) — 

rO)-*'*’ 

4 G v t 

r v J— 

X 

Loading 

Maximum strain 

Error 

Error 

Tension 

0.01 

0.0045 

0.0087 

Tension 

0.10 

0.0273 

0.0719 

Tension 

0.50 

0.1428 

0.4260 

Tension 

1.00 

0.2923 

1.021 

Compression 

-0.01 

-0.0003 

-0.0045 

Compression 

-0.10 

-0.0200 

-0.0591 

Compression 

-0.50 

-0.0334 

-0.1800 
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Table 3. Errors obtained predicting dissipation, Case 2, see Figure 4. 




2 U v (t) 

JA- 

sV 

4 G v t 

CO 

1 

1 II 

Loading 

Maximum strain 

Error 

Error 

Tension 

0.01 

0.0064 

0.0088 

Tension 

0.10 

0.0420 

0.0654 

Tension 

0.50 

0.2450 

0.3890 

Tension 

1.00 

0.5610 

0.9330 

Compression 

-0.01 

-0.0036 

-0.0054 

Compression 

-0.10 

-0.0359 

-0.0567 

Compression 

-0.50 

-0.1110 

-0.1810 


Table 4. Errors obtained predicting dissipation, Case 3, see Figure 5. 






2 U v (t) 

r(t)= 

t 

M S V :S V 
r(t) = 

4 Gy T 

Loading 

Frequency 

Hz 

Minimum 

strain 

Maximum 

strain 

Error 

Error 

Tension 

0.10 

0.004 

0.02 

0.0131 

0.0131 

Tension 

0.10 

0.10 

0.50 

0.3830 

0.3820 

Tension 

1.0 

0.004 

0.02 

0.0132 

0.0132 

Tension 

1.0 

0.10 

0.50 

0.3760 

0.3820 

Tension 

10.0 

0.004 

0.02 

0.0136 

0.0137 

Tension 

10.0 

0.10 

0.50 

0.3720 

0.3840 

Compression 

0.10 

-0.004 

-0.02 

-0.0109 

-0.0109 

Compression 

0.10 

-0.10 

-0.50 

-0.1810 

-0.1830 

Compression 

1.0 

-0.004 

-0.02 

-0.0108 

-0.0108 

Compression 

1.0 

-0.10 

-0.50 

-0.1730 

-0.1830 

Compression 

10.0 

-0.004 

-0.02 

-0.0102 

-0.0104 

Compression 

10.0 

-0.10 

-0.50 

-0.1700 

-0.1840 
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Figure 1. Maxwell solid with one viscous leg. 
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F t 



Plane stress and plane strain 


F t = total force 
u = displacement 


Figure 2. Finite element mesh for uniform axisymmetric, plane stress, and plane strain solids. 
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Total force, F T Displacement, u 



u i 


Displacement, u 


Figure 3. Case 1, load displacement loop. 
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Figure 4. Case 2, load displacement loop. 


t\ - 0.05r 
t 2 - 1\ = 10t 
t 3~ t 2 = 0-05r 

— f4- f 3= 30r 
T = relaxation time constant 
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Total force 









0.7 


Frequency = 0.1 Hz 



Nominal strain 


Figure 8. Dissipation error as a function of nominal strain, frequency = 0.1 Hz. 


Frequency = 0.1 Hz 


0.7 



(I, -3) 


Figure 9. Dissipation error as a function of (/j -3), frequency = 0.1 Hz. 


23 




Axisymmetric FE Model 

Internal heat generated 
by the dissipated 
viscoelastic energy 

Steel disk 
Rubber cylinder 

Heat convects to air 
at outer surface 


Heat conducts to fixture 
at interface 


Dynamic load 

I I I I I I I I I I I 





— 

— 

-1- L 



Eft 


\ 

__ 


3 




H 


Fixture 


Figure 11. Rubber cylinder, finite element mesh, fixture, and steel disk. 
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Load Displacement 



0 


-0.001 


-0.002 

| 

-0.003 

O 

0> 

-0.004 

£ 

<D 

O 

Oh 

-0.005 

-0.006 

s 

-0.007 


-0.008 



Figure 13. Viscoelastic softening of a tall uniform cylinder (H = 0.05 m). 
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Figure 14. Increase in temperature as a function of time at the center of each uniform cylinder. 



Figure 15. Deformation of a tall cylinder with an internal disk (H = 0.05 m). 
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Figure 20. Rubber cylinder with an internal steel ball (H=0.05 m). 


Steel ball 


Boundary conditions are the same as in Figure 13 

▲ 



Figure 21. Deformation of a rubber cylinder with an internal steel ball (FI =0.05 m). 
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Figure 22. Temperature distribution in a tall rubber cylinder with an internal steel ball 
(H=0.05 m). 



Figure 23. Temperature as a function of time for a tall rubber cylinder with an internal steel 
ball (H=0.05 m). 
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